In [1]:
from IPython.core.debugger import Tracer; debug_here = Tracer() #this is the approach that works for ipython debugging
%matplotlib inline
from IPython.display import Image
%config InlineBackend.figure_format = 'svg'

# Export slides with terminal command: 
# jupyter nbconvert run_haats_parallel.ipynb --to slides --post serve --reveal-prefix http://cdn.bootcss.com/reveal.js/3.3.0
# To hide code input cells use
# jupyter nbconvert run_haats_parallel.ipynb --to slides --post serve --reveal-prefix http://cdn.bootcss.com/reveal.js/3.3.0 --template output_toggle
# Export html using  
# jupyter nbconvert run_haats_parallel.ipynb --to html
#then type "ctrl-c" in terminal and then go to html location to re-open slides directly

# Use nbviewer to see results:
# http://nbviewer.jupyter.org/github/GinoAndTonic/ssylvain_public/blob/master/research/haats/run_haats_parallel.ipynb
# http://nbviewer.jupyter.org/github/GinoAndTonic/ssylvain_public/blob/master/research/haats/run_haats_parallel.slides.html    
# or, go to http://nbviewer.jupyter.org/
# Then copy and past the link below into the URL search bar:
# For the notebook:
# https://github.com/GinoAndTonic/ssylvain_public/blob/master/research/haats/run_haats_parallel.ipynb
# For the slides:
# https://github.com/GinoAndTonic/ssylvain_public/blob/master/research/haats/run_haats_parallel.slides.html

from pathos.parallel import ParallelPool as PPool
import copy
import bs4
from import_data import *
from estimation import*
np.random.seed(222)
plt.close("all")
from matplotlib import rc
# import ipdb
import glob

On fancy HA²TS

An introduction to Hidden-state Arbitrage-free Affine Term Structure (HAATS) models

Author: Serginio “Gino” Sylvain

https://github.com/GinoAndTonic/ssylvain_public/blob/master/research/haats https://github.com/GinoAndTonic/ssylvain_public/blob/master/research/haats/haats_documentation.pdf

The US government issues both nominal bonds and inflation-linked bonds (Treasury Inflation Protected Securities; TIPS).

Hence, or a given maturity, the latter tends to be more expensive $\Rightarrow$ lower yields

How do we know when to go long/short Break-Even’s (Nominal Bonds vs ILBs) ?

What is the market-implied inflation rate?

How can we extract the term premium from using a theoretical/structural model?

How can we forecast both nominal and ilb yields in consistent/theoretically robust manner?

How can we extract and forecast market implied discount rates which price these securities?

We often use the Break-Even rate (BE = Nominal Yield - TIPS Yield) as an estimate of expected inflation

For example: if Break-Evens (BE) come down quite a bit, what is that? Should we go long or short Break-evens (by trading ILB’s and short Nominal Bonds or using derivatives)?

But wait…. BE = Exp. Inf. + IRP

  • Perhaps shorting BE’s is more attractive if Exp. Inf. implied by BE’s is “too low”
  • If instead, it is the IRP that is low (or negative) and dragging BE’s down, then it may be reflecting that investors expect future inflation to coincide with a period of higher income growth and/or that nominal yields are coming down sharply due to safe haven flows
  • In general, Exp. Inf. and IRP implied by BE can help inform our investment decisions

Theory

HA²TS ingredients

The model replicates / is heavily inspired by the series of papers by Christensen – Diebold – Lopez – Rudebusch (2007, 2010, 2013)

Suppose there several nominal bonds (N) and several inflation-linked ("real") bonds (R)

The no-arbitrage price of a zero-coupon bond with maturity $\tau$ is:

\begin{eqnarray*} P_{t}^{N}\{\tau\}&=&E_{t}\left[\frac{M_{t+\tau}^{N}}{M_{t}^{N}}\times1\right]=exp\left(-y_{t}^{N}\left\{ \tau\right\} \cdot\tau\right)\\&&\\P_{t}^{R}\{\tau\}&=&E_{t}\left[\frac{M_{t+\tau}^{R}}{M_{t}^{R}}\times1\right]=exp\left(-y_{t}^{R}\left\{ \tau\right\} \cdot\tau\right) \end{eqnarray*}

Here, $y_{t}^{N}$ and $y_{t}^{R}$ are the nominal and real yields respectively. $M_{t}^{N}$ and $M_{t}^{N}$ are the nominal and real state price densities.

The SPDs follow

\begin{eqnarray*} \frac{dM_{t}^{R}}{M_{t}^{R}}&=&-r_{t}^{R}dt-\Gamma_{t}\cdot dW_{t}^{Q}\\&&\\\frac{dM_{t}^{N}}{M_{t}^{N}}&=&-r_{t}^{N}dt-\Gamma_{t}\cdot dW_{t}^{Q} \end{eqnarray*}

For what follows, to simplify the notation, let us suppress the N and R superscripts.

The short rates and risk prices are assumed to be affined in the state variables. This is where we put the rabbit inside the hat...

\begin{eqnarray*} r_{t}&=&\rho_{0}+\rho_{1}\cdot X_{t} \\ \Gamma_{t}&=&\gamma_{0}+\gamma_{1}\cdot X_{t} \end{eqnarray*}
\begin{eqnarray*} P_{t}\{\tau\}&=&E_{t}^{Q}\left[exp\left(-\int_{t}^{t+\tau}r_{s}ds\right)\right]\\&&\\r_{t}&=&\rho_{0}+\rho_{1}\cdot X_{t} \\&&\\\Rightarrow P_{t}\{\tau\}&=&exp\left(B_{t}\left\{ \tau\right\} \cdot X_{t}+G_{t}\left\{ \tau\right\} \right)\\&&\text{Since }exp\left(-\int_{0}^{t}r_{s}ds\right)P_{t}\left\{ \tau\right\} \text{ is Martingale under }Q\text{, it's drift is zero.}\\&&\text{Thus, we allow for no arbitrage and }\\&&\text{ by using Ito's Lemman and setting }E_{t}^{Q}\left[\frac{dP_{t}\{\tau\}}{P_{t}\{\tau\}}\right]-r_{t}=0\\&&B_{t}\left\{ \tau\right\} \text{ and }G_{t}\left\{ \tau\right\} \text{ solve some ODEs.}\\&&\\y_{t}\{\tau\}&=&-\frac{1}{\tau}ln\left(P_{t}\{\tau\}\right)=-\frac{1}{\tau}B_{t}\left\{ \tau\right\} \cdot X_{t}-\frac{1}{\tau}G_{t}\left\{ \tau\right\} \end{eqnarray*}

Solving for B{t} and G{t} uniquely typically requires imposing some ad-hoc parameter values and other restrictions with little motivation.

Instead, following Christensen, Lopez, Rudebusch (2010) we assume a dynamic Nelson-Seigel (1987) model and impose level, slope and curvature restrictions by replacing the Nelson-Seigel (1987) parameters with the level, slope and curvature state variables.

This is sensible since there is a lot of evidence that Level, Slope, and Curvature (e.g. from PCA) explain the cross-section of government bonds. Furthermore, the data we will use to fit the model will itself be smoothe yields data from Nelson-Seigel-Svennson-type models.

Hence, the second key assumption concerns the state variables:

\begin{eqnarray*} X_{t}&=&\left(\begin{array}{c} L_{t}^{N}\\ S_{t}\\ C_{t}\\ L_{t}^{R} \end{array}\right) \\ dX_{t}&=&K^{P}\left(\theta^{P}-X_{t}\right)dt+\Sigma dW_{t} \end{eqnarray*}

This implies

\begin{eqnarray*} y_{t}^{N}\{\tau\} &=& L_{t}^{N}+S_{t}\left(\frac{1-e^{-\lambda\tau}}{\lambda\tau}\right)+C_{t}\left(\frac{1-e^{-\lambda\tau}}{\lambda\tau}-e^{-\lambda\tau}\right)-\frac{G_{t}^{N}\left\{ \tau\right\} }{\tau} \\ y_{t}^{R}\{\tau\} &=& L_{t}^{R}+\alpha^{R}S_{t}\left(\frac{1-e^{-\lambda\tau}}{\lambda\tau}\right)+\alpha^{R}C_{t}\left(\frac{1-e^{-\lambda\tau}}{\lambda\tau}-e^{-\lambda\tau}\right)-\frac{G_{t}^{R}\left\{ \tau\right\} }{\tau} \end{eqnarray*}

Empirical approach

We observe the Nominal bond and ILB yields but the state variables (X) are hidden. We also need to estimate the model parameters.

First, we can re-write the key equations of the model...

Measurements: $$y_{t}=A_{0}+A_{1}X_{t}+\epsilon_{t} \qquad \text{with }\epsilon_{t}\sim N(0,\Phi)$$

States: $$X_{t}=U_0+U_1 X_{t_-1}+\eta_{t} \qquad \text{with } \eta_{t}\sim N(0,Q)$$

Since the Brownian Motion increments are Gaussian, Kalman-filtering is an efficient and consistent estimator. It also allows for asynchronous variables (crucial: if we later include observed Macro variables as state variables).

We use the Kalman filter along with Expectation-Maximization (EM) algorithm to jointly extract the hidden states and estimate the model parameters.

Note that $y_t$ is both the input and the target.

EM Algorithm:

1) Start with a guess for the set of parameters, $\Omega$

2) Run the Kalman filter

3) Run the Kalman smoother

4) Compute the expected log-likelihood:

$$ \mathcal{\tilde{L}}\left\{ \Omega\right\} = E\left[ln\left(\left(\prod_{t=1}^{T}p\left(Y_{t}|X_{t},\Omega\right)\right)\left(\prod_{t=1}^{T}p\left(X_{t}|X_{t-1},\Omega\right)\right)\right)\bigg|X^{(T)}\right] $$

5) Solve for $\hat{\Omega}$ $$\hat{\Omega} = arg\max_{\Omega}\mathcal{\tilde{L}}\left\{ \Omega\right\} $$

6) repeat steps 2-5 until change in $\hat{\Omega}$ is below a pre-specified tolerance level

Bayesian EM Algorithm: replace steps 4. and 5. with

4) choosing priors for $\Omega$ and

5) sampling from posterior distribution for $\Omega$

Data

Importing (Nelson Siegel smoothed/fitted) yield data for TIPS and Nominal bonds...

Data source:

http://www.federalreserve.gov/econresdata/researchdata/feds200628.xls

https://www.federalreserve.gov/econresdata/researchdata/feds200805.xls

In [2]:
# Import data from FED
tips_data, nominal_data = ImportData.importUS_Data(plots=1,save=1)
processing time for importing US ILB data: 16.5218708334
            TIPSY02  TIPSY03  TIPSY04  TIPSY05  TIPSY06  TIPSY07  TIPSY08  \
1999-01-04      NaN      NaN      NaN   3.9244   3.9369   3.9310   3.9170   
1999-01-05      NaN      NaN      NaN   3.9385   3.9525   3.9475   3.9340   
1999-01-06      NaN      NaN      NaN   3.8522   3.8357   3.8182   3.8026   
1999-01-07      NaN      NaN      NaN   3.8995   3.8827   3.8550   3.8282   
1999-01-08      NaN      NaN      NaN   3.9043   3.8985   3.8771   3.8525   

            TIPSY09  TIPSY10  TIPSY11  TIPSY12  TIPSY13  TIPSY14  TIPSY15  \
1999-01-04   3.9002   3.8832   3.8673   3.8527   3.8397   3.8281   3.8179   
1999-01-05   3.9175   3.9005   3.8845   3.8698   3.8566   3.8449   3.8345   
1999-01-06   3.7893   3.7783   3.7691   3.7613   3.7547   3.7490   3.7441   
1999-01-07   3.8051   3.7859   3.7698   3.7562   3.7448   3.7349   3.7264   
1999-01-08   3.8293   3.8089   3.7913   3.7762   3.7633   3.7521   3.7424   

            TIPSY16  TIPSY17  TIPSY18  TIPSY19  TIPSY20  
1999-01-04   3.8088   3.8007   3.7934   3.7869   3.7810  
1999-01-05   3.8252   3.8169   3.8096   3.8029   3.7969  
1999-01-06   3.7398   3.7360   3.7326   3.7296   3.7269  
1999-01-07   3.7189   3.7123   3.7064   3.7012   3.6965  
1999-01-08   3.7338   3.7263   3.7196   3.7136   3.7082  
           TIPSY02      TIPSY03      TIPSY04      TIPSY05      TIPSY06  \
count  3160.000000  3160.000000  3160.000000  4407.000000  4407.000000   
mean      0.104398     0.245663     0.411247     1.247476     1.388270   
std       1.449942     1.329237     1.249470     1.571417     1.514691   
min      -2.202000    -2.067000    -1.892800    -1.724000    -1.538900   
25%            NaN          NaN          NaN          NaN          NaN   
50%            NaN          NaN          NaN          NaN          NaN   
75%            NaN          NaN          NaN          NaN          NaN   
max       5.479800     4.486200     4.013800     4.378400     4.385600   

           TIPSY07      TIPSY08      TIPSY09      TIPSY10      TIPSY11  \
count  4407.000000  4407.000000  4407.000000  4407.000000  4407.000000   
mean      1.507833     1.610125     1.697962     1.773473     1.838364   
std       1.457489     1.402522     1.351277     1.304448     1.262255   
min      -1.357800    -1.188300    -1.016600    -0.862700    -0.724600   
25%            NaN          NaN          NaN          NaN          NaN   
50%            NaN          NaN          NaN          NaN          NaN   
75%            NaN          NaN          NaN          NaN          NaN   
max       4.383800     4.380900     4.377400     4.371100     4.366800   

           TIPSY12      TIPSY13      TIPSY14      TIPSY15      TIPSY16  \
count  4407.000000  4407.000000  4407.000000  4407.000000  4407.000000   
mean      1.894053     1.941741     1.982471     2.017155     2.046591   
std       1.224629     1.191346     1.162083     1.136486     1.114183   
min      -0.598300    -0.484300    -0.382600    -0.292300    -0.212700   
25%            NaN          NaN          NaN          NaN          NaN   
50%            NaN          NaN          NaN          NaN          NaN   
75%            NaN          NaN          NaN          NaN          NaN   
max       4.363600     4.360900     4.358500     4.356500     4.354700   

           TIPSY17      TIPSY18      TIPSY19      TIPSY20  
count  4407.000000  4407.000000  4407.000000  4407.000000  
mean      2.071483     2.092451     2.110042     2.124740  
std       1.094812     1.078029     1.063517     1.050984  
min      -0.142600    -0.080900    -0.026800     0.020800  
25%            NaN          NaN          NaN          NaN  
50%            NaN          NaN          NaN          NaN  
75%            NaN          NaN          NaN          NaN  
max       4.353100     4.351700     4.350500     4.349400  
tips_data column kept are:
Index([u'TIPSY02', u'TIPSY03', u'TIPSY04', u'TIPSY05', u'TIPSY06', u'TIPSY07',
       u'TIPSY08', u'TIPSY09', u'TIPSY10', u'TIPSY11', u'TIPSY12', u'TIPSY13',
       u'TIPSY14', u'TIPSY15', u'TIPSY16', u'TIPSY17', u'TIPSY18', u'TIPSY19',
       u'TIPSY20'],
      dtype='object')
c:\anaconda2\lib\site-packages\numpy\lib\function_base.py:3403: RuntimeWarning: Invalid value encountered in median
  RuntimeWarning)
processing time for importing US NB data: 39.1022234726
            SVENY01  SVENY02  SVENY03  SVENY04  SVENY05  SVENY06  SVENY07  \
1961-06-14   2.9825   3.3771   3.5530   3.6439   3.6987   3.7351   3.7612   
1961-06-15   2.9941   3.4137   3.5981   3.6930   3.7501   3.7882   3.8154   
1961-06-16   3.0012   3.4142   3.5994   3.6953   3.7531   3.7917   3.8192   
1961-06-19   2.9949   3.4386   3.6252   3.7199   3.7768   3.8147   3.8418   
1961-06-20   2.9833   3.4101   3.5986   3.6952   3.7533   3.7921   3.8198   

            SVENY08  SVENY09  SVENY10   ...     SVENY21  SVENY22  SVENY23  \
1961-06-14      NaN      NaN      NaN   ...         NaN      NaN      NaN   
1961-06-15      NaN      NaN      NaN   ...         NaN      NaN      NaN   
1961-06-16      NaN      NaN      NaN   ...         NaN      NaN      NaN   
1961-06-19      NaN      NaN      NaN   ...         NaN      NaN      NaN   
1961-06-20      NaN      NaN      NaN   ...         NaN      NaN      NaN   

            SVENY24  SVENY25  SVENY26  SVENY27  SVENY28  SVENY29  SVENY30  
1961-06-14      NaN      NaN      NaN      NaN      NaN      NaN      NaN  
1961-06-15      NaN      NaN      NaN      NaN      NaN      NaN      NaN  
1961-06-16      NaN      NaN      NaN      NaN      NaN      NaN      NaN  
1961-06-19      NaN      NaN      NaN      NaN      NaN      NaN      NaN  
1961-06-20      NaN      NaN      NaN      NaN      NaN      NaN      NaN  

[5 rows x 30 columns]
            SVENY01       SVENY02       SVENY03       SVENY04       SVENY05  \
count  13739.000000  13739.000000  13739.000000  13739.000000  13739.000000   
mean       5.221512      5.450177      5.633840      5.787301      5.918971   
std        3.281879      3.200316      3.108261      3.019517      2.938882   
min        0.082800      0.158300      0.302500      0.431300      0.588600   
25%             NaN           NaN           NaN           NaN           NaN   
50%             NaN           NaN           NaN           NaN           NaN   
75%             NaN           NaN           NaN           NaN           NaN   
max       16.462000     15.911800     15.574600     15.349800     15.177600   

            SVENY06       SVENY07       SVENY08       SVENY09       SVENY10  \
count  13739.000000  13739.000000  11203.000000  11203.000000  11203.000000   
mean       6.033364      6.133223      6.478916      6.570095      6.649857   
std        2.867758      2.806193      2.937031      2.879881      2.831142   
min        0.762000      0.941900      1.121400      1.295500      1.400700   
25%             NaN           NaN           NaN           NaN           NaN   
50%             NaN           NaN           NaN           NaN           NaN   
75%             NaN           NaN           NaN           NaN           NaN   
max       15.061100     15.010900     14.979900     14.957200     14.939800   

          ...           SVENY21      SVENY22      SVENY23      SVENY24  \
count     ...       7649.000000  7649.000000  7649.000000  7649.000000   
mean      ...          5.941984     5.947573     5.949433     5.948006   
std       ...          1.919230     1.911577     1.904812     1.898851   
min       ...          1.896900     1.945100     1.995000     2.046600   
25%       ...               NaN          NaN          NaN          NaN   
50%       ...               NaN          NaN          NaN          NaN   
75%       ...               NaN          NaN          NaN          NaN   
max       ...         10.485600    10.520600    10.555300    10.589300   

          SVENY25      SVENY26      SVENY27      SVENY28      SVENY29  \
count  7649.00000  7649.000000  7649.000000  7649.000000  7649.000000   
mean      5.94369     5.936842     5.927777     5.916780     5.904104   
std       1.89362     1.889061     1.885132     1.881798     1.879028   
min       2.09970     2.154300     2.210300     2.267400     2.325700   
25%           NaN          NaN          NaN          NaN          NaN   
50%           NaN          NaN          NaN          NaN          NaN   
75%           NaN          NaN          NaN          NaN          NaN   
max      10.62170    10.652600    10.681900    10.709900    10.736500   

           SVENY30  
count  7649.000000  
mean      5.889974  
std       1.876806  
min       2.351800  
25%            NaN  
50%            NaN  
75%            NaN  
max      10.761800  

[8 rows x 30 columns]
nominal_data column kept are:
Index([u'SVENY01', u'SVENY02', u'SVENY03', u'SVENY04', u'SVENY05', u'SVENY06',
       u'SVENY07', u'SVENY08', u'SVENY09', u'SVENY10', u'SVENY11', u'SVENY12',
       u'SVENY13', u'SVENY14', u'SVENY15', u'SVENY16', u'SVENY17', u'SVENY18',
       u'SVENY19', u'SVENY20', u'SVENY21', u'SVENY22', u'SVENY23', u'SVENY24',
       u'SVENY25', u'SVENY26', u'SVENY27', u'SVENY28', u'SVENY29', u'SVENY30'],
      dtype='object')
In [3]:
# The data is daily but let's focus of sub-sample of 5yrs with weekly frequency
tips_data, nominal_data = tips_data.loc['2004-01-01':'2016-05-30'],nominal_data.loc['2004-01-01':'2016-05-30']
tips_data, nominal_data = tips_data.asfreq('W', method='pad'), nominal_data.asfreq('W', method='pad')
In [4]:
#import seaborn as sns
import seaborn.apionly as sns #use sns.distplot but maintain the default matplotlib styling
sns.set("talk", font_scale=1, rc={"lines.linewidth": 1,"axes.labelcolor":'black',"text.color":'black'})
In [5]:
fig = plt.figure();ax1= plt.subplot(121)
nominal_data.plot(ax=ax1)
plt.legend(loc='center left',fontsize=9,frameon=0, bbox_to_anchor=(1, 0.5))

ax2=plt.subplot(122,sharey=ax1)
tips_data.plot(ax=ax2)
plt.legend(loc='center left',fontsize=9,frameon=0, bbox_to_anchor=(1, 0.5))
fig.subplots_adjust(wspace=0.5);plt.show()
In [6]:
nominal_data[['Nominals_y01','Nominals_y05','Nominals_y10','Nominals_y30']].dropna(0).describe()
Out[6]:
Nominals_y01 Nominals_y05 Nominals_y10 Nominals_y30
count 646.000000 646.000000 646.000000 646.000000
mean 1.528356 2.492177 3.443919 4.082252
std 1.739018 1.298076 1.078136 0.782272
min 0.091700 0.617300 1.478600 2.351800
25% 0.216300 1.498075 2.378675 3.320150
50% 0.465700 2.096400 3.637900 4.287750
75% 2.678700 3.659175 4.408500 4.697850
max 5.277600 5.110100 5.256600 5.824800
In [7]:
tips_data[['TIPS_y02','TIPS_y05','TIPS_y10','TIPS_y20']].dropna(0).describe()
Out[7]:
TIPS_y02 TIPS_y05 TIPS_y10 TIPS_y20
count 647.000000 647.000000 647.000000 647.000000
mean 0.119568 0.589371 1.194858 1.637228
std 1.460355 1.192668 0.957266 0.727053
min -2.202000 -1.712700 -0.830400 0.071500
25% -1.049900 -0.326300 0.495650 1.000800
50% -0.128200 0.506500 1.358800 1.900000
75% 0.976600 1.440500 1.959000 2.245800
max 5.159100 3.679600 3.577300 3.268200

There is a fair amout of correlation and auto-correlation.

The model can handle this (does not blow up). But, we may want to be more careful about which bonds we choose to include

In [8]:
plt.rc('text', usetex=False);fig = plt.figure(figsize=(20,8))
ax1=plt.subplot(121)
sns.heatmap(    nominal_data.corr())

ax2=plt.subplot(122);plt.rc('text', usetex=False)
sns.heatmap(    tips_data.corr())
fig.subplots_adjust(wspace=0.3);plt.show()
In [9]:
from itertools import cycle
fig = plt.figure(figsize=(16,6))

ax1=plt.subplot(121);lines = ["-","--","-.",":"];linecycler = cycle(lines)
for i in nominal_data.columns:
    plt.acorr(nominal_data[[i]].dropna(0).iloc[:,0],maxlags=50, usevlines=False, linestyle=next(linecycler),marker=None,label=i )
    plt.xlim(0,50)
plt.title('Nominal Bonds autocorrelations');plt.xlabel('lag')
plt.legend(loc='center left',fontsize=9,frameon=0, bbox_to_anchor=(1, 0.5))

ax2=plt.subplot(122,sharey=ax1);lines = ["-","--","-.",":"];linecycler = cycle(lines)
for i in tips_data.columns:
    plt.acorr(tips_data[[i]].dropna(0).iloc[:,0],maxlags=50, usevlines=False, linestyle=next(linecycler),marker=None,label=i )
    plt.xlim(0,50)
plt.title('TIPS autocorrelations');plt.xlabel('lag')
plt.legend(loc='center left',fontsize=9,frameon=0, bbox_to_anchor=(1, 0.5))
fig.subplots_adjust(wspace=0.5);plt.show()

Running code...

In [10]:
start_time = time.time()

# The interveal between each rolling window: the gap by which the estimationd window shifts
# (e.g. with tgap = 1, rolling window is updated daily)
tgap = 30

# Rolling window: 0 if using expanding window, 1 if using rolling window
rolling = 0

#Rolling window size: size of rolling window (in years).. Use inf for full sample estimation
windowsize = np.inf;

np.set_printoptions(precision=32, suppress=True) #increase precision on  numeric values



################################################

# PRIMITIVES:
figures = []
# use allow_missing_data= 1 to extract ILB and Nominal dates where both are non-missing
allow_missing_data = 0

# set frequency of the data: daily, monthly, quarterly, yearly
estim_freq = 'weekly'

fix_Phi = 0     # "1" if you want to fix the volatility of observed yields using covar of historical data
                # "0" if you want to jointly estimate it with other model parameters
setdiag_Kp = 0  # "1" if you want to Kp to be diagonal so the state variables are assumed independent
                # "0" if you want to Kp to be unrestricted

# options for initializing the Kalman filter error variance:
#'steady_state' or 'unconditional' or 'identity' matrix
initV = 'unconditional'

# number of hidden state variables 4, or 6
num_states = 4

# Specify the maturities of data we want to use
US_ilbmaturities = np.array([2, 3,  5, 6, 8, 9, 10])
US_nominalmaturities = np.array([2, 3,  5, 6, 8, 9, 10])
US_maturities = np.hstack((US_nominalmaturities, US_ilbmaturities))

############################################################

# Set start and end dates for estimation
sdate, edate = '2004-01-01', '2016-05-30'#time.strftime("%Y-%m-%d") #'2010-01-01'
print("start date: %s" % sdate)
print("end date: %s" % edate)

# extract data for desired maturities and dates
tips_data, nominal_data = ImportData.importUS_Data(US_ilbmaturities, US_nominalmaturities,plots=0,save=1)
data = ImportData.extract_subset(tips_data, nominal_data, sdate, edate, allow_missing_data, estim_freq)
start date: 2004-01-01
end date: 2016-05-30
In [11]:
# Instantiate estimation object:
estimation1 =Rolling()
In [12]:
# Set up data, etc.:
estimation1.run_setup(data, US_ilbmaturities, US_nominalmaturities, \
                estim_freq=estim_freq, num_states=num_states,\
                fix_Phi=fix_Phi, setdiag_Kp=setdiag_Kp, initV=initV)
In [13]:
# estimation_method,tolerance, maxiter ,toltype,solver_mle,maxiter_mle, maxfev_mle, ftol_mle, xtol_mle, \
#     constraints_mle, priors_bayesian, maxiter_bayesian, burnin_bayesian, multistart = 'em_mle', 1e-4, 2 , 'max_abs', \
#     'Nelder-Mead', 5, 5, 0.0001, 0.0001, 'off', None, 1000, None, 4,8

estimation_method,tolerance, maxiter ,toltype,solver_mle,maxiter_mle, maxfev_mle, ftol_mle, xtol_mle, \
    constraints_mle, priors_bayesian, maxiter_bayesian, burnin_bayesian, multistart, ncpus = 'em_mle', 1e-6, 100 , \
    'max_abs', 'Powell', 10000, 10000, 0.01, 0.01, 'off', None, 1000, None, 100, 8

estimation1.fit(estimation_method=estimation_method, tolerance=tolerance, maxiter=maxiter, toltype=toltype, \
                solver_mle=solver_mle, maxiter_mle=maxiter_mle, maxfev_mle=maxfev_mle, ftol_mle=ftol_mle,
                xtol_mle=xtol_mle, constraints_mle=constraints_mle, \
                priors_bayesian=priors_bayesian, maxiter_bayesian=maxiter_bayesian, burnin_bayesian=burnin_bayesian,
                multistart=multistart, ncpus=ncpus)
done running parallel workers
Cannot kill worker with PID: 13680
('Produced error:', <class 'psutil.NoSuchProcess'>)
Cannot kill worker with PID: 12560
('Produced error:', <class 'psutil.NoSuchProcess'>)
Cannot kill worker with PID: 8328
('Produced error:', <class 'psutil.NoSuchProcess'>)
Cannot kill worker with PID: 8252
('Produced error:', <class 'psutil.NoSuchProcess'>)
Cannot kill worker with PID: 10628
('Produced error:', <class 'psutil.NoSuchProcess'>)
Cannot kill worker with PID: 9344
('Produced error:', <class 'psutil.NoSuchProcess'>)
Cannot kill worker with PID: 8200
('Produced error:', <class 'psutil.NoSuchProcess'>)
Cannot kill worker with PID: 14844
('Produced error:', <class 'psutil.NoSuchProcess'>)
Job execution statistics:
 job count | % of all jobs | job time sum | time per job | job server
       108 |        100.00 |     3542039.8830 |  32796.665583 | local
Time elapsed since server creation 470730.459
0 active tasks, 8 cores


Out[13]:
(array([ 0.7189983057298259261358452931745 ,
         0.06675021655694526878654926349554,
        -0.05477296333168651942813198729709,
        -0.660791830183465545012211350695  ,
         2.56458769755862014960712258471176,
         0.31136181496603315688176394360198,
        -0.02415388135657227722252748947085,
        -2.64621339178163772487550886580721,
        -1.89648535947129626144658232078655,
         0.41258793039867480700877422350459,
         2.24240542706412337992105676676147,
         4.12060167435922064527176189585589,
        -1.50481260346075007561239544884302,
        -0.19140504311887512889356344203406,
         0.01505359510811383333406254791953,
         1.78391281527102063542145060637267,
         0.0000000000000136006162496472636 ,
         0.00000014695850269930723083503725,
         0.00000002673012205457334445580422,
         0.00000000000001343137217161515797,
         0.00000001218438236693340827846711,
         0.00000000000001344129538628759155,
         0.00000004865898711082913095599345,
         0.00002451655708907563550489824844,
         0.00000871057689471139874463120195,
         0.00000060997124987088177500274338,
         0.00000000000001343665734582748412,
         0.00000155131939898648308883148846,
         0.00000297414353120477706991374618,
         0.00000452321855531813098974698584,
         0.76903071148298129955378499289509,
         0.51686018623532947735554898827104,
         0.01293542781880345840517243516388,
         0.0244559864121003320280856030422 ,
         0.05729202727317606047563458560035,
         0.0122999716909977518941587959489 ,
         0.09630303848304933023172225148301,
        -0.22193345808191020096522549920337,
        -0.03824844683127729544347417345307,
         0.03580641254058725825304421164219]),
  final_simplex: (array([[ 0.71827016915181363110320944542764,
         0.06837780152128108968057063066226,
        -0.05477560930172503672608286251489, ...,
        -0.22196701425774431615245418925042,
        -0.03825532309271722358756662174528,
         0.03581176177052522446864202265715],
       [ 0.7182701691232629137573439948028 ,
         0.06837780152175693126892497275549,
        -0.05477560930141131545489940890548, ...,
        -0.22196701425950832176070548484859,
        -0.03825532309256370749883657822465,
         0.03581176177080383493667170569097],
       [ 0.71827016904235296834713153657503,
         0.0683778015228913849110625733374 ,
        -0.05477560930051311727240204163536, ...,
        -0.22196701424526776857959475819371,
        -0.03825532309245961715138406589176,
         0.03581176176842684744094924553792],
       ..., 
       [ 0.7182701684229488847677203011699 ,
         0.06837780152242722842004241101677,
        -0.05477560930189823151792438693519, ...,
        -0.22196701430055401638163914412871,
        -0.03825532309572508143924451928797,
         0.03581176177113992720180135620467],
       [ 0.71827016807998500524945484357886,
         0.06837780154018074929567205799685,
        -0.05477560929186701954574090223105, ...,
        -0.22196701425855994149749506050284,
        -0.03825532310632848154163809795136,
         0.03581176176621050921689004553627],
       [ 0.71827016783492503293473419034854,
         0.06837780149016643715942365133742,
        -0.05477560929433698821977571924435, ...,
        -0.22196701428253240440113813747303,
        -0.03825532308820012855621328640154,
         0.03581176176064758409633981273146]]), array([-86205.27050047021475620567798614501953,
       -86205.27050046654767356812953948974609,
       -86205.27050045937357936054468154907227,
       -86205.27050045572104863822460174560547,
       -86205.27050045420764945447444915771484,
       -86205.27050045337819028645753860473633,
       -86205.27050045299984049052000045776367,
       -86205.27050044981297105550765991210938,
       -86205.27050044466159306466579437255859,
       -86205.27050043993222061544656753540039,
       -86205.27050043552299030125141143798828,
       -86205.27050043530471157282590866088867,
       -86205.27050043475173879414796829223633,
       -86205.27050043235067278146743774414062,
       -86205.27050043235067278146743774414062,
       -86205.27050043175404425710439682006836,
       -86205.27050042808696161955595016479492,
       -86205.27050042776681948453187942504883,
       -86205.2705004257877590134739875793457 ,
       -86205.27050042556948028504848480224609,
       -86205.27050042452174238860607147216797,
       -86205.27050042056362144649028778076172,
       -86205.27050041370966937392950057983398,
       -86205.27050041258917190134525299072266,
       -86205.2705004107410786673426628112793 ,
       -86205.27050040860194712877273559570312,
       -86205.27050040647736750543117523193359,
       -86205.2705004048330010846257209777832 ,
       -86205.27050040051108226180076599121094,
       -86205.27050039827008731663227081298828,
       -86205.27050039492314681410789489746094,
       -86205.27050039230380207300186157226562,
       -86205.27050038971356116235256195068359,
       -86205.27050038540619425475597381591797,
       -86205.27050038131710607558488845825195,
       -86205.27050037418666761368513107299805,
       -86205.27050036544096656143665313720703,
       -86205.27050035487627610564231872558594,
       -86205.27050034439889714121818542480469,
       -86205.27050028793746605515480041503906,
       -86205.27049988575163297355175018310547]))
           fun: -86212.586536269635
       message: 'Optimization terminated successfully.'
          nfev: 2993
           nit: 1608
        status: 0
       success: True
             x: array([  0.7189983057298259261358452931745 ,
         0.06675021655694526878654926349554,
        -0.05477296333168651942813198729709,
        -0.660791830183465545012211350695  ,
         2.56458769755862014960712258471176,
         0.31136181496603315688176394360198,
        -0.02415388135657227722252748947085,
        -2.64621339178163772487550886580721,
        -1.89648535947129626144658232078655,
         0.41258793039867480700877422350459,
         2.24240542706412337992105676676147,
         4.12060167435922064527176189585589,
        -1.50481260346075007561239544884302,
        -0.19140504311887512889356344203406,
         0.01505359510811383333406254791953,
         1.78391281527102063542145060637267,
       -31.92866129072119463216949952766299,
       -15.73311558458282100048108986811712,
       -17.43747474060423741093472926877439,
       -31.94118321747151156841937336139381,
       -18.22311083913640672449218982364982,
       -31.9404446815519129643234919058159 ,
       -16.83842931545032328699562640395015,
       -10.61616186910917569719003950012848,
       -11.65097253568685076174915593583137,
       -14.30985401224898190264411823591217,
       -31.94078980013682667049579322338104,
       -13.37640476396721567198255797848105,
       -12.72555444908524258096349512925372,
       -12.306286747760871946866245707497  ,
         0.76903071148298129955378499289509,
        -0.65998287386093124062824699649354,
        -4.34778538937981018364098417805508,
        -3.71088024999973065121139370603487,
        -2.85959380502651017152970780443866,
        -4.39815831815131197402024554321542,
         0.09630303848304933023172225148301,
        -0.22193345808191020096522549920337,
        -0.03824844683127729544347417345307,
         0.03580641254058725825304421164219]),
 <estimation.Rolling instance at 0x00000000136F7308>)
In [14]:
estimation1.fit_path.tail()
Out[14]:
objective criteria Kp[0] Kp[1] Kp[2] Kp[3] Kp[4] Kp[5] Kp[6] Kp[7] ... a[0] lmda[0] sigmas[0] sigmas[1] sigmas[2] sigmas[3] thetap[0] thetap[1] thetap[2] thetap[3]
iter
11 -86171.893116 0.118138 0.740713 0.066369 -0.054187 -0.662338 2.703451 0.317151 -0.020856 -2.756052 ... 0.769031 0.51686 0.012935 0.024456 0.057292 0.0123 0.095733 -0.220797 -0.035352 0.036163
12 -86196.214511 0.107629 0.737272 0.066635 -0.054469 -0.661105 2.647902 0.315157 -0.021231 -2.710904 ... 0.769031 0.51686 0.012935 0.024456 0.057292 0.0123 0.096031 -0.221132 -0.036290 0.036034
13 -86205.488369 0.100085 0.728073 0.066650 -0.054485 -0.661071 2.601204 0.313239 -0.022432 -2.673853 ... 0.769031 0.51686 0.012935 0.024456 0.057292 0.0123 0.096251 -0.221464 -0.037234 0.035914
14 -86212.589473 0.096311 0.718998 0.066750 -0.054773 -0.660792 2.564588 0.311362 -0.024154 -2.646213 ... 0.769031 0.51686 0.012935 0.024456 0.057292 0.0123 0.096303 -0.221933 -0.038248 0.035806
15 -86212.586536 0.000000 0.718998 0.066750 -0.054773 -0.660792 2.564588 0.311362 -0.024154 -2.646213 ... 0.769031 0.51686 0.012935 0.024456 0.057292 0.0123 0.096303 -0.221933 -0.038248 0.035806

5 rows × 42 columns

In [15]:
# Delete temporary files:
map(os.remove, glob.glob(r""+str.replace(os.getcwd(), '\\', '/')+"/output/parallel_worker_output"+"*.txt"))
Out[15]:
[None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None,
 None]
In [16]:
fig = plt.figure(figsize=(14,8));ax1=plt.subplot(221);estimation1.fit_path['objective'].plot(ax=ax1)
plt.title('optimality: negative log-likelihood');plt.xlabel('outer iterations')

ax2 = plt.subplot(222);estimation1.fit_path_inner['sub_objective'].plot(ax=ax2)
plt.title('optimality: negative log-likelihood');plt.xlabel('inner iterations')
plt.ylim(ymax=estimation1.fit_path_inner['sub_objective'].iloc[1],ymin=estimation1.fit_path_inner['sub_objective'].min())

ax3 = plt.subplot(212);estimation1.fit_path['criteria'].plot()
plt.title('tolerance: maximum absolute change in parameters');fig.subplots_adjust(hspace=0.3);plt.show();

Examining some results

Model fit

We do a decent job at fitting some bond yields but a not so good job at fitting others.

This can be improved upon by having more flexible model (more parameters) and a more robust optimization method.

In [17]:
estimation1.collect_results()
processing time for forecast: 8.36822667054
In [18]:
linestyles = ['-', '--', '-.', ':','-', '--', '-.', ':','-', '--', '-.', ':','-', '--', '-.', ':']
plt.rc('text', usetex=True);fig = plt.figure(figsize=(16,8));ax1 = plt.subplot(121)
pd.concat([estimation1.Y.iloc[3:,0:4].rename(columns={i:str.replace(i,'_','\_') for i in estimation1.Y.columns.values},inplace=False),
    estimation1.Ytt_new.iloc[3:,0:4].rename(columns={i:str.replace(i,'_','\_') +'$_{k|k}$' for i in estimation1.Ytt_new.columns.values},inplace=False),
    estimation1.Yttl_new.iloc[3:,0:4].rename(columns={i:str.replace(i,'_','\_') +'$_{k|k-1}$' for i in estimation1.Yttl_new.columns.values},inplace=False)]
, axis=1).plot(ax=ax1,style=linestyles,linewidth=1)
plt.legend(loc='best',fontsize=9,frameon=0);plt.title('Measurements and Filtered Measurements');plt.axes.labelcolor='black'

linestyles = ['-', '--', '-.','-', '--', '-.', '-', '--', '-.','-', '--', '-.']
ax2 = plt.subplot(122)
pd.concat([estimation1.Y.iloc[3:,4:7].rename(columns={i:str.replace(i,'_','\_') for i in estimation1.Y.columns.values},inplace=False),
    estimation1.Ytt_new.iloc[3:,4:7].rename(columns={i:str.replace(i,'_','\_') +'$_{k|k}$' for i in estimation1.Ytt_new.columns.values},inplace=False),
    estimation1.Yttl_new.iloc[3:,4:7].rename(columns={i:str.replace(i,'_','\_') +'$_{k|k-1}$' for i in estimation1.Yttl_new.columns.values},inplace=False)]
, axis=1).plot(ax=ax2,style=linestyles,linewidth=1)
plt.legend(loc='best',fontsize=9,frameon=0);plt.title('Measurements and Filtered Measurements');plt.axes.labelcolor='black';plt.show();
In [19]:
linestyles = ['-', '--', '-.', ':','-', '--', '-.', ':','-', '--', '-.', ':','-', '--', '-.', ':']
plt.rc('text', usetex=True);fig = plt.figure(figsize=(16,8));ax1 = plt.subplot(121)
pd.concat(
    [estimation1.Y.iloc[3:,7:11].rename(columns={i:str.replace(i,'_','\_') for i in estimation1.Y.columns.values},inplace=False),
    estimation1.Ytt_new.iloc[3:,7:11].rename(columns={i:str.replace(i,'_','\_') +'$_{k|k}$' for i in estimation1.Ytt_new.columns.values},inplace=False),
    estimation1.Yttl_new.iloc[3:,7:11].rename(columns={i:str.replace(i,'_','\_') +'$_{k|k-1}$' for i in estimation1.Yttl_new.columns.values},inplace=False)]
, axis=1).plot(ax=ax1,style=linestyles,linewidth=1)
plt.legend(loc='best',fontsize=9,frameon=0);plt.title('Measurements and Filtered Measurements');plt.axes.labelcolor='black'

linestyles = ['-', '--', '-.','-', '--', '-.','-', '--', '-.','-', '--', '-.'];ax2 = plt.subplot(122)
pd.concat(
    [estimation1.Y.iloc[3:,11:].rename(columns={i:str.replace(i,'_','\_') for i in estimation1.Y.columns.values},inplace=False),
    estimation1.Ytt_new.iloc[3:,11:].rename(columns={i:str.replace(i,'_','\_') +'$_{k|k}$' for i in estimation1.Ytt_new.columns.values},inplace=False),
    estimation1.Yttl_new.iloc[3:,11:].rename(columns={i:str.replace(i,'_','\_') +'$_{k|k-1}$' for i in estimation1.Yttl_new.columns.values},inplace=False)]
, axis=1).plot(ax=ax2,style=linestyles,linewidth=1)
plt.legend(loc='best',fontsize=9,frameon=0);plt.title('Measurements and Filtered Measurements')
plt.axes.labelcolor='black';plt.show()
In [20]:
fit_e, fit_se, fit_mse_all, fit_rmse_all =  estimation1.fit_e, \
                        estimation1.fit_se, estimation1.fit_mse_all, \
                        estimation1.fit_rmse_all 
In [21]:
plt.rc('text', usetex=False);fig = plt.figure(figsize=(14,7));ax1 = plt.subplot(221)
fit_se.iloc[:,0:7].plot(ax=ax1,linewidth=1 )
plt.legend(loc='center left',fontsize=9,frameon=0, bbox_to_anchor=(1, 0.5))
plt.title('Squared-Error for Nominal Bonds at each date',fontsize=12)
plt.axes.labelcolor='black'

ax2 = plt.subplot(222);fit_se.iloc[:,7:].plot(ax=ax2,linewidth=1 )
plt.legend(loc='center left',fontsize=9,frameon=0, bbox_to_anchor=(1, 0.5))
plt.title('Squared-Error for TIPS at each date',fontsize=12);plt.axes.labelcolor='black'

ax3 = plt.subplot(212)
fit_rmse_all.plot(ax=ax3,linewidth=1,kind='bar')
plt.legend(loc='best',fontsize=9,frameon=0)
plt.title('RMSE across all dates',fontsize=12)
plt.axes.labelcolor='black';fig.subplots_adjust(wspace=0.5,hspace=0.3);plt.show() 

Implied States/Latent Variables

The first values for the filtered states is usually off... the value from the smoother are not...

In [22]:
linestyles = ['-', '--', '-.', ':','-', '--', '-.', ':','-', '--', '-.', ':','-', '--', '-.', ':'];plt.rc('text', usetex=True)
fig, ax = plt.subplots(1)
figures = {'fig2': fig, 'ax_fig2': ax}
pd.concat(
    [estimation1.XtT_new.rename(columns={i:i+'$_{smoother}$' for i in estimation1.Xtt_new.columns.values},inplace=False),
    estimation1.Xtt_new.rename(columns={i:i+'$_{k|k}$' for i in estimation1.Xtt_new.columns.values},inplace=False),
    estimation1.Xttl_new.rename(columns={i:i+'$_{k|k-1}$' for i in estimation1.Xttl_new.columns.values},inplace=False)]
, axis=1).plot(ax=figures['ax_fig2'],figsize=(8,8),style=linestyles,linewidth=1)
plt.legend(loc='center left',fontsize=9,frameon=0, bbox_to_anchor=(1, 0.5))
figures['ax_fig2'].set_title('Filtered States/Latent Variables');plt.axes.labelcolor='black';plt.show()

Implied Expected Inflation and Probability of Deflation

We are not yet confident in the implied inflation and deflation probabilities because the fit is not yet satisfactory...

Nontheless, to illustrate the machinery we plot the results below.

In [23]:
estimation1.expected_inflation()
In [26]:
plt.rc('text', usetex=False);fig, ax = plt.subplots(1)
figures = {'fig2': fig, 'ax_fig2': ax}
estimation1.exp_inf[['horizon_1yr','horizon_2yr','horizon_5yr','horizon_10yr','horizon_30yr']].plot(ax=figures['ax_fig2'],
                                                                                                    figsize=(6,6),linewidth=1)
plt.legend(loc='center left',fontsize=12,frameon=0, bbox_to_anchor=(1, 0.5))
figures['ax_fig2'].set_title('Implied Expected Inflation')
plt.axes.labelcolor='black';plt.show()
In [27]:
plt.rc('text', usetex=False);fig, ax = plt.subplots(1);figures = {'fig2': fig, 'ax_fig2': ax}
estimation1.prob_def[['horizon_1yr','horizon_2yr','horizon_5yr','horizon_10yr','horizon_30yr']].plot(ax=figures['ax_fig2'],
                                                                                                    figsize=(6,6),linewidth=1)
plt.legend(loc='center left',frameon=0, bbox_to_anchor=(1, 0.5))
figures['ax_fig2'].set_title('Implied Probability of Deflation')
plt.axes.labelcolor='black';plt.show()

Forecasting (in-sample)

In the forecast below, we are “cheating” a bit because the parameters and filtered states are obtained for a full-sample fit.

Instead what we should do is to use a rolling or expanding window to do the estimation and then forecast out of sample starting from the last date of the estimation window.

Nonetheless, for now let us examine some in-sample forecasts

Not surprisingly, these forecasts are rather poor given that the fit itself was poor...

In [28]:
forecast_e, forecast_se, forecast_mse, forecast_rmse, forecast_mse_all, forecast_rmse_all =  estimation1.forecast_e, \
                        estimation1.forecast_se, estimation1.forecast_mse, estimation1.forecast_rmse, estimation1.forecast_mse_all, \
                        estimation1.forecast_rmse_all
        
forecast = estimation1.yields_forecast 
forecast_std = estimation1.yields_forecast_std       
In [29]:
fig, ax = sns.plt.subplots(1,figsize=(14,8))
line,=sns.plt.plot(forecast[['Nominals_y05']].iloc[np.arange(0, forecast.shape[0], estimation1.forecast_horizon+1)].set_index(\
    forecast[['Nominals_y05']].iloc[np.arange(0, forecast.shape[0], estimation1.forecast_horizon+1)].index.get_level_values('date')),linewidth=1.5\
                  )
line.set_label('Nominals_y05: actual')
for t in forecast.index.get_level_values('date').unique():#loop and plot forecast for each horizon
    line,=sns.plt.plot(
        forecast[['Nominals_y05']].iloc[forecast.index.get_level_values('date')==t,:].set_index(
        forecast[['Nominals_y05']].iloc[forecast.index.get_level_values('date')==t,:].index.get_level_values('horizon'))
        ,color='red',linewidth=0.5)
    fill=plt.fill_between(forecast[['Nominals_y05']].iloc[forecast.index.get_level_values('date')==t,:].index.get_level_values('horizon').values
                     , 
                     forecast[['Nominals_y05']].iloc[forecast.index.get_level_values('date')==t,:].values[:,0]-
                     forecast_std[['Nominals_y05']].iloc[forecast.index.get_level_values('date')==t,:].values[:,0],
                     forecast[['Nominals_y05']].iloc[forecast.index.get_level_values('date')==t,:].values[:,0]+
                     forecast_std[['Nominals_y05']].iloc[forecast.index.get_level_values('date')==t,:].values[:,0]
                    , facecolor='red', interpolate=True, alpha=.05
                    )
line.set_label('Nominals_y05: forecasts');fill.set_label('Nominals_y05: forecasts stderr')
sns.plt.legend(loc='best',fontsize=9);sns.plt.axes.labelcolor='black';sns.plt.show()
In [30]:
fig, ax = sns.plt.subplots(1,figsize=(14,8))
line,=sns.plt.plot(forecast[['TIPS_y08']].iloc[np.arange(0, forecast.shape[0], estimation1.forecast_horizon+1)].set_index(\
    forecast[['TIPS_y08']].iloc[np.arange(0, forecast.shape[0], estimation1.forecast_horizon+1)].index.get_level_values('date')),linewidth=1.5\
                  )
line.set_label('TIPS_y08: actual')
for t in forecast.index.get_level_values('date').unique():#loop and plot forecast for each horizon
    line,=sns.plt.plot(
        forecast[['TIPS_y08']].iloc[forecast.index.get_level_values('date')==t,:].set_index(
        forecast[['TIPS_y08']].iloc[forecast.index.get_level_values('date')==t,:].index.get_level_values('horizon'))
        ,color='red',linewidth=0.5)
    fill=plt.fill_between(forecast[['TIPS_y08']].iloc[forecast.index.get_level_values('date')==t,:].index.get_level_values('horizon').values
                     , 
                     forecast[['TIPS_y08']].iloc[forecast.index.get_level_values('date')==t,:].values[:,0]-
                     forecast_std[['TIPS_y08']].iloc[forecast.index.get_level_values('date')==t,:].values[:,0],
                     forecast[['TIPS_y08']].iloc[forecast.index.get_level_values('date')==t,:].values[:,0]+
                     forecast_std[['TIPS_y08']].iloc[forecast.index.get_level_values('date')==t,:].values[:,0]
                    , facecolor='red', interpolate=True, alpha=.05
                    )
line.set_label('TIPS_y08: forecasts');fill.set_label('TIPS_y08: forecasts stderr')
sns.plt.legend(loc='best',fontsize=9);sns.plt.axes.labelcolor='black';sns.plt.show()
In [31]:
plt.rc('text', usetex=False);fig = plt.figure(figsize=(14,7));ax1 = plt.subplot(221)
forecast_rmse.iloc[:,0:7].plot(ax=ax1,linewidth=1)
plt.legend(loc='center left',fontsize=9,frameon=0, bbox_to_anchor=(1, 0.5))
plt.title('RMSE of forecasts for Nominal Bonds at each date',fontsize=12);plt.axes.labelcolor='black'

ax2 = plt.subplot(222);forecast_rmse.iloc[:,7:].plot(ax=ax2,linewidth=1 )
plt.legend(loc='center left',fontsize=9,frameon=0, bbox_to_anchor=(1, 0.5))
plt.title('RMSE of forecasts for TIPS at each date',fontsize=12);plt.axes.labelcolor='black'

ax3 = plt.subplot(212);forecast_rmse_all.plot(ax=ax3,linewidth=1,kind='bar' )
plt.legend(loc='best',fontsize=9,frameon=0);plt.title('RMSE of forecasts across all dates',fontsize=12);plt.axes.labelcolor='black';
fig.subplots_adjust(wspace=0.5,hspace=0.3);plt.show()
In [ ]:
from IPython.core.debugger import Tracer; debug_here = Tracer() #this is the approach that works for ipython debugging

import estimation 
reload(estimation)
from estimation import *

import kalman 
reload(kalman)
from kalman import *



# import estim_constraints 
# reload(estim_constraints)
# from estim_constraints import *
In [32]:
import copy
estimation2 = copy.deepcopy(estimation1)

Conclusion

Take away

  • This is a decent start...
  • There is a lot more to be done...

Lingering questions

  • Why not use a simpler approach (like PCA) to extract the latent state variables?
    • Answer: because no arbritrage condition would be violated by such a model and if we would want to use such a model to identify arbitrage, we would be unable to do so because it would not be clear whether we ideed identify arbitrage opportunities of if the model's flaws drive the results. Furthermore, we would not be able to forecast future yields.
  • Why not use VAR?
    • Answer: again, no arbritrage condition would be violated